**CALCULATES UTILITY FROM DIFFERENT MAJOR CHOICES VIA SIMULATION OF A CRRA UTILITY MODEL
set more off
set seed 984575
pause on
				
*SET GLOBALS FOR KEY PARAMETERS


	*SET WORKING TIME TO MATCH AGE 22 THROUGH 37 (15 YEARS/60 quarters - number is 59 b/c initial period is t = 0)
	global working_time = 79 
	global numquarters = $working_time - 20 + 1
	
	*NUMBER OF SIMULATION RUNS
	global run = 200
	



***FOUR - YEARS***	
	
	*GROWTH RATES FOR EACH FIELD (FROM ESTIMATES)
		local growth_lib = 38.663
		local growth_ag = `growth_lib' + 11.961
		local growth_com = `growth_lib' + 3.001
		local growth_it = `growth_lib' +  26.781
		local growth_voc = `growth_lib' +  12.256
		local growth_eng = `growth_lib' +  35.451
		local growth_bio = `growth_lib' +  40.746
		local growth_sci = `growth_lib' +  21.979
		local growth_soc = `growth_lib' +   6.329
		local growth_bus = `growth_lib' +  28.267
		local growth_und = `growth_lib' +  14.000
	
	*SD (FROM ESTIMATES)
		local sd_lib = 3603.68
		local sd_ag = `sd_lib' + 668.449
		local sd_com = `sd_lib' + 454.274
		local sd_it = `sd_lib' + 1063.187
		local sd_voc = `sd_lib' + 774.428
		local sd_eng = `sd_lib' + 2251.211
		local sd_bio = `sd_lib' + 1693.051
		local sd_sci = `sd_lib' + 1002.743
		local sd_soc = `sd_lib' + 350.938
		local sd_bus = `sd_lib' + 1731.414
		local sd_und = `sd_lib' + 804.240
				

	*STARTING VALUES FOR EACH FIELD (FROM "ALPHA" ESTIMATES UNDER MEAN) 
		local alpha_lib = 5564.723
		local alpha_ag = `alpha_lib'  + 124.126
		local alpha_com = `alpha_lib' + 317.128
		local alpha_it = `alpha_lib'  + 2154.616
		local alpha_voc = `alpha_lib' + 952.730
		local alpha_eng = `alpha_lib' + 3506.859
		local alpha_bio = `alpha_lib' + 120.482
		local alpha_sci = `alpha_lib' + 433.165
		local alpha_soc = `alpha_lib' + -495.106
		local alpha_bus = `alpha_lib' + 2573.439
		local alpha_und = `alpha_lib' + -2200.990
	
				
*LOAD QTE DATA
use "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\qte_bootstrap_data_4yr.dta", clear
cap postclose simulation_4yr
postfile simulation_4yr str40(major sector) float(discountrate) int(run) float(gamma) int(quantiles decilenum) str20(scenario) float(utility_major utility_libarts relutility)  using "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_4yr.dta" , replace

		

*EVERYTHING IS RELATIVE TO LIBERAL ARTS
keep if relmajor == "lib_4yr"


*LOOP OVER DISCOUNT RATE

	
*DISCOUNT RATE (QUARTERLY)
foreach beta of numlist .9975 .995 .99 {

	*INTEREST RATE (QUARTERLY)
	global r = 1-`beta'
	
*LOOP OVER COEFFICIENT OF RELATIVE RISK AVERSION
 foreach gamma of numlist  0 0.75 1.5 2 3{
 di ""
 di "gamma = " `gamma'


   *LOOP THROUGH MAJORS
   foreach major in "ag" "bio" "bus" "com" "eng" "it" "sci" "soc" "und" "voc"  {		
		di ""
		di "Major = `major'"
		
		*LOOP OVER DRAWS OF SIMULATION
		forvalues run=1/$run {
			di "Run =" `run'

	   
		*CONVERT WAGE ESTIMATES TO MATRIX FORM
		mkmat quantile major_wage relmajor_wage if major == "`major'_4yr", matrix(quantiles)
		
		*SET PROBABILITY OF LANDING IN THE QUANTILE DISTRIBUTION
			local percentile = runiform(0.5,9.5)
			forvalues decile = 1(1)9 {
				if `percentile' >= `decile'-.5 & `percentile'<`decile'+.5 local major_wage = quantiles[`decile',2]
				if `percentile' >= `decile'-.5 & `percentile'<`decile'+.5 local relmajor_wage = quantiles[`decile',3]
				if `percentile' >= `decile'-.5 & `percentile'<`decile'+.5 local decilenum = `decile'
			}
			*di "major_wage = " `major_wage'
			*di "relmajor_wage = " `relmajor_wage'
		
		*CALCULATE STARTING VALUE FOR ASSIGNED QUANTILE - REMEMBER, START AT PERIOD 20 WHICH IS 5 YEARS AFTER HS GRADUATION
			local beta_sum_major = 0
			local beta_sum_relmajor = 0
				forvalues t = 20/$working_time {
				  local beta_sum_major = `beta_sum_major'+`growth_`major''*`t'	
				  local beta_sum_relmajor = `beta_sum_relmajor'+`growth_lib' *`t'	
				}
			local beta_sum_major = `beta_sum_major'/($working_time - 20 + 1)
			local beta_sum_relmajor = `beta_sum_relmajor'/($working_time - 20 + 1)
			local alpha_major = `major_wage' - `beta_sum_major'
			local alpha_relmajor = `relmajor_wage' - `beta_sum_relmajor' 
			
			*di "quantile = " `percentile'
			*di "major_wage = " `major_wage'
			*di "relmajor_wage = " `relmajor_wage'
			*di "beta_sum_major = " `beta_sum_major'
			*di "beta_sum_relmajor = " `beta_sum_relmajor'
			*di "alpha_major = " `alpha_major'
			*di "alpha_relmajor = " `alpha_relmajor'
			
			

	
		
		matrix input temp = (.,.,.,.,.,.,.,.,.,.,.)
		matrix colnames temp = "quarter" "mean major - expected" "mean major - realized" "quantile major - expected"  ///
			"quantile major - realized"  "" "mean relmajor - expected" "mean relmajor - realized" /// 
			"quantiles relmajor - expected" "quantiles relmajor - realized" ""
	
		*LOOP OVER QUARTERS STARTING AT 20
		forvalues t=20/$working_time {

			
			*CALCULATE PREDICTED EARNINGS IN EACH PERIOD - REMEMBER, START AT PERIOD 20 WHICH IS 5 YEARS AFTER HS GRADUATION
			*DRAW EARNINGS FROM DISTRIBUTION OF MEAN = PREDICTED ERANINGS AND ESTIMATED SD FOR THE MAJOR
			*SET NEGATIVE VALUES TO ZERO AS UTILITY MODEL REQUIES POSITIVE EARNINGS UNDER SOME VALUES OF GAMMA

				*WITH QUANTILES
				local major_earnings =  max(rnormal(`alpha_major' + `growth_`major''*`t', `sd_`major''),0)
				local major_earnings_expected = `alpha_major' + `growth_`major''*`t'
				local relmajor_earnings =  max(rnormal(`alpha_relmajor' + `growth_lib'*`t',`sd_lib'),0)
				local relmajor_earnings_expected =  `alpha_relmajor' + `growth_lib'*`t'

				*MEAN EARNINGS ONLY
							
				local major_earnings_mean =  max(rnormal(`alpha_`major'' + `growth_`major''*`t', `sd_`major''),0)
				local major_earnings_expected_mean = `alpha_`major'' + `growth_`major''*`t'
				local relmajor_earnings_mean =  max(rnormal(`alpha_lib' + `growth_lib'*`t',`sd_lib'),0)
				local relmajor_earnings_expected_mean =  `alpha_lib' + `growth_lib'*`t'

					
				*di "major_earnings_mean = " `major_earnings_mean'
				*di "relmajor_earnings_mean = " `relmajor_earnings_mean'
		
				
			matrix  temp = temp\ (`t', `major_earnings_expected_mean', `major_earnings_mean', `major_earnings_expected', `major_earnings',. , ///
				`relmajor_earnings_expected_mean', `relmajor_earnings_mean', `relmajor_earnings_expected', `relmajor_earnings',. )
			}
			
			
			
			matrix predicted_earnings = temp[2...,1...]
		
		
			

*DO EVERYTHING WITH QUANTILES AND W/O
		
		
		
	*CALCULATE UTILITY OF FIELD AND REFERENCE FIELD
		
		*UNCERTAINTY WITH NO SAVINGS, Mean
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',3] >0 local utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',3]^(1-`gamma'))/(1-`gamma')
			if predicted_earnings[`t',8] >0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',8]^(1-`gamma') )/(1-`gamma')
		 }
		local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')
		post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (0) (.) ("Uncertainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')

				
		 *UNCERTAINTY WITH NO SAVINGS, QUANTILES
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',5] >0  local utility_major = `utility_major'+ `beta'^(`t')*((predicted_earnings[`t',5]^(1-`gamma') )/(1-`gamma'))
			if predicted_earnings[`t',10] >0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*((predicted_earnings[`t',10]^(1-`gamma') )/(1-`gamma'))
		 }
		local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')
		  post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (1) (`decilenum') ("Uncertainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
				
				
		*CERTAINTY WITH NO SAVINGS (NO CV,  INCOME STREAM IS DETERMINITIVE), MEAN
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',2] >0  local utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',2]^(1-`gamma') )/(1-`gamma')
			if predicted_earnings[`t',7] >0  local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',7]^(1-`gamma') )/(1-`gamma')
			
		 }
		local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')
		  post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (0) (.) ("Certainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
		 
		 
		*CERTAINTY WITH NO SAVINGS (NO CV,  INCOME STREAM IS DETERMINITIVE), QUANTILES
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',4] >0  local utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',4]^(1-`gamma') )/(1-`gamma')
			if predicted_earnings[`t',9] >0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',9]^(1-`gamma') )/(1-`gamma')
		 }
		local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')
		 post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (1) (`decilenum') ("Certainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
		
	
		
		*PERMANENT INCOME HYPOTHESIS - INCOME HAS SAME PDV (DISCOUNTED TO t=20 IN EACH PERIOD
		
				  *CALCULATE PDV OF EARNINGS
						
					*CREATE MATRIX TO HOLD CONSUMPTION VALUES UNDER PIH
					matrix input pih = (.,.,.,.,.)
					matrix colnames pih = "quarter" "mean major" "quantile major" "mean relmajor" "quantile relmajor" 
						
				
						*MAJOR, MEAN
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',3]/((1 + $r)^`t')
					   }
					   local pdv_major_expected_mean_avg = `pdv_major_expected'/$numquarters
					   
					  *RELMAJOR, MEAN
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',8]/((1 + $r)^`t')
					   }
					   local pdv_relmajor_expected_mean_avg = `pdv_major_expected'/$numquarters
					   
					   *MAJOR, QUANTILES
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',5]/((1 + $r)^`t')
					   }
					   local pdv_major_expected_quant_avg = `pdv_major_expected'/$numquarters
					   
					   
					   *RELMAJOR, QUANTILES
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',10]/((1 + $r)^`t')
					   }
					   local pdv_relmajor_expected_quant_avg = `pdv_major_expected'/$numquarters
					 
					   
					   *PLACE VALUES OF EQUAL PDV IN EACH QUARTER IN A MATRIX
					   forvalues t = 1/$numquarters{
						  	local pih_mean_major = `pdv_major_expected_mean_avg'/((1-$r)^`t')
						  	local pih_mean_relmajor = `pdv_relmajor_expected_mean_avg'/((1-$r)^`t')
						  	local pih_quant_major = `pdv_major_expected_quant_avg'/((1-$r)^`t')
						  	local pih_quant_relmajor = `pdv_relmajor_expected_quant_avg'/((1-$r)^`t')
							
							
							matrix  pih = pih\(`t' + 20, `pih_mean_major',`pih_quant_major',`pih_mean_relmajor',`pih_quant_relmajor')
					   }

						matrix pih = pih[2...,1...]
				   
				   
				   
				   *CALCULATE UTILITY - MEAN
					 local utility_major = 0
					 local utility_relmajor = 0

					 forvalues t = 1/$numquarters {
						if pih[`t',2] > 0  local utility_major = `utility_major'+ `beta'^(`t')*(pih[`t',2]^(1-`gamma') )/(1-`gamma')
						if pih[`t',4] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(pih[`t',4]^(1-`gamma') )/(1-`gamma')		
					 }
					local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')	
					post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (0) (.) ("Permanent Income Hypothesis") (`utility_major') (`utility_relmajor') (`rel_utility')

					 
				   *CALCULATE UTILITY - QUANTILES
					 local utility_major = 0
					 local utility_relmajor = 0

					 forvalues t = 1/$numquarters {
					    if pih[`t',3] > 0 	local utility_major = `utility_major'+ `beta'^(`t')*(pih[`t',3]^(1-`gamma') )/(1-`gamma')
						if pih[`t',5] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(pih[`t',5]^(1-`gamma') )/(1-`gamma')
					 }
		local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')
					post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (1) (`decilenum') ("Permanent Income Hypothesis") (`utility_major') (`utility_relmajor') (`rel_utility')		   
			
 
				   
		*SAVINGS RULE - SAVE ANYTHING ABOVE KNOWN INCOME STREAM, DRAW DOWN SAVINGS UNTIL HIT PREDICTED INCOME OR SAVINGS = 0
		
			 *MEAN
			 local utility_major = 0
			 local utility_relmajor = 0
			 local savings_major = 0
			 local savings_relmajor = 0
			 forvalues t = 1/$numquarters {
				if `savings_major' >= 0 { 
					 local consumption_major = min(predicted_earnings[`t',2], predicted_earnings[`t',3] + `savings_major')
				}
				else {
					 local consumption_major = predicted_earnings[`t',3]
				}
			    local savings_major = (predicted_earnings[`t',3] + `savings_major' - `consumption_major')*(1+$r) 
				
				if `savings_relmajor' >= 0 { 
					 local consumption_relmajor = min(predicted_earnings[`t',7], predicted_earnings[`t',8] + `savings_major')
				}
				else {
					 local consumption_relmajor = predicted_earnings[`t',8]
				}
			    local savings_relmajor = (predicted_earnings[`t',8] + `savings_relmajor' - `consumption_relmajor')*(1+$r) 

				
				if `consumption_major' > 0 local utility_major = `utility_major'+ `beta'^(`t')*(`consumption_major'^(1-`gamma') )/(1-`gamma')
				if `consumption_relmajor' > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(`consumption_relmajor'^(1-`gamma') )/(1-`gamma')
				
				/*
				di `t'
				di `savings_major'
				di `consumption_major'
				di `utility_major'
				di ""
				di `savings_relmajor'
				di `consumption_relmajor'
				di `utility_relmajor'
				*/
			}
				
			 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')	
			 *di `rel_utility'
			 post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (0) (.) ("Savings Rule") (`utility_major') (`utility_relmajor') (`rel_utility')
			
			 
			
			 *QUANTILE
			 local utility_major = 0
			 local utility_relmajor = 0
			 local savings_major = 0
			 local savings_relmajor = 0
			 
			  forvalues t = 1/$numquarters {
				if `savings_major' >= 0 { 
					 local consumption_major = min(predicted_earnings[`t',4], predicted_earnings[`t',5] + `savings_major')
				}
				else {
					 local consumption_major = predicted_earnings[`t',5]
				}
			    local savings_major = (predicted_earnings[`t',5] + `savings_major' - `consumption_major')*(1+$r) 
				
				if `savings_relmajor' >= 0 { 
					 local consumption_relmajor = min(predicted_earnings[`t',9], predicted_earnings[`t',10] + `savings_major')
				}
				else {
					 local consumption_relmajor = predicted_earnings[`t',10]
				}
			    local savings_relmajor = (predicted_earnings[`t',10] + `savings_relmajor' - `consumption_relmajor')*(1+$r) 
	
				if `consumption_major' > 0 local utility_major = `utility_major'+ `beta'^(`t')*(`consumption_major'^(1-`gamma') )/(1-`gamma')
				if `consumption_relmajor' > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(`consumption_relmajor'^(1-`gamma') )/(1-`gamma')
			}
				
			 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
			 post simulation_4yr ("`major'") ("4yr") (`beta') (`run') (`gamma') (1) (`decilenum') ("Savings Rule") (`utility_major') (`utility_relmajor') (`rel_utility')

	
	*CLOSE LOOP FOR RUN
	}
	

*CLOSE LOOP FOR MAJOR
}

*CLOSE LOOP FOR CRRA (GAMMA)
} 		

*CLOSE LOOP FOR dISCOUNT RATE (BETA)
}
		postclose simulation_4yr
		use "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_4yr.dta" , clear
		save "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_4yr_runs.dta" , replace
	

		*COLLAPSE TO SUMMARY DATA BY MAJOR AND TYPE
		collapse (median) utility_major utility_libarts relutility, by(major scenario gamma quantiles discountrate)
		sort major scenario quantiles  gamma 
		matrix drop _all
		save "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_4yr.dta" , replace

*/
		

***TWO - YEARS***	
	
	*GROWTH RATES FOR EACH FIELD (FROM ESTIMATES)
		local growth_lib = 29.767
		local growth_ag = `growth_lib' + 3.193
		local growth_com = `growth_lib' + -5.985
		local growth_it = `growth_lib' +  -3.192
		local growth_voc = `growth_lib' +  -1.243
		local growth_eng = `growth_lib' +  -5.944
		local growth_bio = `growth_lib' +  1.753
		local growth_sci = `growth_lib' +  7.475
		local growth_soc = `growth_lib' +  -3.537
		local growth_bus = `growth_lib' +  0.915
		local growth_educ = `growth_lib' + -5.647
		local growth_und = `growth_lib' +  6.065

	
	*SD (FROM ESTIMATES)
		local sd_lib = 3019.146
		local sd_ag = `sd_lib' + 312.055
		local sd_com = `sd_lib' - 255.553
		local sd_it = `sd_lib' - 162.216
		local sd_voc = `sd_lib' + 359.028
		local sd_eng = `sd_lib' + 377.265
		local sd_bio = `sd_lib' + 522.740
		local sd_sci = `sd_lib' + 260.319
		local sd_soc = `sd_lib' - 141.234
		local sd_bus = `sd_lib' + 224.596
		local sd_educ = `sd_lib' - 212.913
		local sd_und = `sd_lib' + 330.116

				

	*STARTING VALUES FOR EACH FIELD (FROM "ALPHA" ESTIMATES UNDER MEAN) 
		local alpha_lib = 3701.744

		local alpha_ag = `alpha_lib'  + 473.634
		local alpha_com = `alpha_lib' + -310.973
		local alpha_it = `alpha_lib'  + 42.687
		local alpha_voc = `alpha_lib' + 1399.865
		local alpha_eng = `alpha_lib' + 747.775
		local alpha_bio = `alpha_lib' + 1054.388
		local alpha_sci = `alpha_lib' + -330.596
		local alpha_soc = `alpha_lib' + -99.185
		local alpha_bus = `alpha_lib' + 682.156
		local alpha_educ = `alpha_lib' + 83.193
		local alpha_und = `alpha_lib' + 256.395


				
*LOAD QTE DATA
use "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\qte_bootstrap_data_2yr.dta", clear
cap postclose simulation_2yr
postfile simulation_2yr str40(major sector) float(discountrate) int(run) float(gamma) int(quantiles) str20(scenario) double(utility_major utility_libarts relutility)  using "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_2yr.dta" , replace

*EVERYTHING IS RELATIVE TO LIBERAL ARTS
keep if relmajor == "lib_2yr"

*DISCOUNT RATE (QUARTERLY)
foreach beta of numlist .9975 .995 .99 {

	*INTEREST RATE (QUARTERLY)
	global r = 1-`beta'

*LOOP OVER COEFFICIENT OF RELATIVE RISK AVERSION
 foreach gamma of numlist  0 0.75 1.5 2 3 {
 di ""
 di "gamma = " `gamma'

	
	

   *LOOP THROUGH MAJORS
   foreach major in "ag" "bio" "bus" "com" "educ" "eng" "it" "sci" "soc" "und" "voc" {
		di ""
		di "Major = `major'"
		
		*LOOP OVER DRAWS OF SIMULATION
		forvalues run=1/$run {
			di "Run =" `run'

	   
		*CONVERT WAGE ESTIMATES TO MATRIX FORM
		mkmat quantile major_wage relmajor_wage if major == "`major'_2yr", matrix(quantiles)
		
		*SET PROBABILITY OF LANDING IN THE QUANTILE DISTRIBUTION
			local percentile = runiform(0.5,9.5)
			forvalues decile = 1(1)9 {
				if `percentile' >= `decile'-.5 & `percentile'<`decile'+.5 local major_wage = quantiles[`decile',2]
				if `percentile' >= `decile'-.5 & `percentile'<`decile'+.5 local relmajor_wage = quantiles[`decile',3]
			}
		
		*CALCULATE STARTING VALUE FOR ASSIGNED QUANTILE - REMEMBER, START AT PERIOD 20 WHICH IS 5 YEARS AFTER HS GRADUATION
			local beta_sum_major = 0
			local beta_sum_relmajor = 0
				forvalues t = 20/$working_time {
				  local beta_sum_major = `beta_sum_major'+`growth_`major''*`t'	
				  local beta_sum_relmajor = `beta_sum_relmajor'+`growth_lib' *`t'	
				}
			local beta_sum_major = `beta_sum_major'/($working_time - 20 + 1)
			local beta_sum_relmajor = `beta_sum_relmajor'/($working_time - 20 + 1)
			local alpha_major = `major_wage' - `beta_sum_major'
			local alpha_relmajor = `relmajor_wage' - `beta_sum_relmajor' 


	
		
		matrix input temp = (.,.,.,.,.,.,.,.,.,.,.)
		matrix colnames temp = "quarter" "mean major - expected" "mean major - realized" "quantile major - expected"  ///
			"quantile major - realized"  ""  "mean relmajor - expected" "mean relmajor - realized" /// 
			"quantiles relmajor - expected" "quantiles relmajor - realized" ""
	
		*LOOP OVER QUARTERS STARTING AT 20
		forvalues t=20/$working_time {


			*CALCULATE PREDICTED EARNINGS IN EACH PERIOD - REMEMBER, START AT PERIOD 20 WHICH IS 5 YEARS AFTER HS GRADUATION
			*DRAW EARNINGS FROM DISTRIBUTION OF MEAN = PREDICTED ERANINGS AND ESTIMATED SD FOR THE MAJOR
			*SET NEGATIVE VALUES TO ZERO AS UTILITY MODEL REQUIES POSITIVE EARNINGS UNDER SOME VALUES OF GAMMA

				
				
				*WITH QUANTILES
				local major_earnings =  max(rnormal(`alpha_major' + `growth_`major''*`t', `sd_`major''),0)
				local major_earnings_expected = `alpha_major' + `growth_`major''*`t'
				local relmajor_earnings =  max(rnormal(`alpha_relmajor' + `growth_lib'*`t',`sd_lib'),0)
				local relmajor_earnings_expected =  `alpha_relmajor' + `growth_lib'*`t'

				*MEAN EARNINGS ONLY
							
				local major_earnings_mean =  max(rnormal(`alpha_`major'' + `growth_`major''*`t', `sd_`major''),0)
				local major_earnings_expected_mean = `alpha_`major'' + `growth_`major''*`t'
				local relmajor_earnings_mean =  max(rnormal(`alpha_lib' + `growth_lib'*`t',`sd_lib'),0)
				local relmajor_earnings_expected_mean =  `alpha_lib' + `growth_lib'*`t'
				
				
			matrix  temp = temp\ (`t', `major_earnings_expected_mean', `major_earnings_mean', `major_earnings_expected', `major_earnings', ., ///
				`relmajor_earnings_expected_mean', `relmajor_earnings_mean', `relmajor_earnings_expected', `relmajor_earnings', .)
			}
			
			
			
			matrix predicted_earnings = temp[2...,1...]
		
		

*DO EVERYTHING WITH QUANTILES AND W/O
		*NOTE I FOUND A STATA MACHINE PRECISINO ISSUE WHERE IT WILL NOT TAKE ROOTS OF ZEROS SO HAVE TO DO AN IF-STATEMENT WORK AROUND, SAI 1-31-24
		
		
	*CALCULATE UTILITY OF FIELD AND REFERENCE FIELD
		
		*UNCERTAINTY WITH NO SAVINGS, Mean
		 local utility_major = 0
 		 local utility_relmajor = 0		 

		 forvalues t = 1/$numquarters {
		 	if predicted_earnings[`t',3] > 0 local  utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',3]^(1-`gamma'))/(1-`gamma')
			if predicted_earnings[`t',8] > 0 local  utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',8]^(1-`gamma') )/(1-`gamma')
		 }
		 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
		 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (0) ("Uncertainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
		

		if `utility_major' == . {
			return
			di `utility_major'
			matrix list predicted_earnings
		}
		
		 *UNCERTAINTY WITH NO SAVINGS, QUANTILES
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',5] > 0  local utility_major = `utility_major'+ `beta'^(`t')*((predicted_earnings[`t',5]^(1-`gamma') )/(1-`gamma'))
			if predicted_earnings[`t',10] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*((predicted_earnings[`t',10]^(1-`gamma') )/(1-`gamma'))
		 }
		 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
		 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (1) ("Uncertainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
		
	
						
		*CERTAINTY WITH NO SAVINGS (NO CV,  INCOME STREAM IS DETERMINITIVE), MEAN
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',2] > 0 local utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',2]^(1-`gamma') )/(1-`gamma')
			if predicted_earnings[`t',7] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',7]^(1-`gamma') )/(1-`gamma')
		 }
		 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
		 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (0) ("Certainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
	
		
		
		*CERTAINTY WITH NO SAVINGS (NO CV,  INCOME STREAM IS DETERMINITIVE), QUANTILES
		 local utility_major = 0
 		 local utility_relmajor = 0

		 forvalues t = 1/$numquarters {
			if predicted_earnings[`t',4] > 0  local utility_major = `utility_major'+ `beta'^(`t')*(predicted_earnings[`t',4]^(1-`gamma') )/(1-`gamma')
			if predicted_earnings[`t',9] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(predicted_earnings[`t',9]^(1-`gamma') )/(1-`gamma')
		 }
		 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
		 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (1) ("Certainty, No Savings") (`utility_major') (`utility_relmajor') (`rel_utility')
	
		
		
		
		*PERMANENT INCOME HYPOTHESIS - INCOME HAS SAME PDV (DISCOUNTED TO t=20 IN EACH PERIOD
		
				  *CALCULATE PDV OF EARNINGS
						
					*CREATE MATRIX TO HOLD CONSUMPTION VALUES UNDER PIH
					matrix input pih = (.,.,.,.,.)
					matrix colnames pih = "quarter" "mean major" "quantile major" "mean relmajor" "quantile relmajor" 
						
				
						*MAJOR, MEAN
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',3]/((1 + $r)^`t')
					   }
					   local pdv_major_expected_mean_avg = `pdv_major_expected'/$numquarters
					   
					  *RELMAJOR, MEAN
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',8]/((1 + $r)^`t')
					   }
					   local pdv_relmajor_expected_mean_avg = `pdv_major_expected'/$numquarters
					   
					   *MAJOR, QUANTILES
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',5]/((1 + $r)^`t')
					   }
					   local pdv_major_expected_quant_avg = `pdv_major_expected'/$numquarters
					   
					   
					   *RELMAJOR, QUANTILES
					   local pdv_major_expected = 0
					   forvalues t = 1/$numquarters{
						 local pdv_major_expected = `pdv_major_expected' + predicted_earnings[`t',10]/((1 + $r)^`t')
					   }
					   local pdv_relmajor_expected_quant_avg = `pdv_major_expected'/$numquarters
					 
					   
					   *PLACE VALUES OF EQUAL PDV IN EACH QUARTER IN A MATRIX
					   forvalues t = 1/$numquarters{
						  	local pih_mean_major = `pdv_major_expected_mean_avg'/((1-$r)^`t')
						  	local pih_mean_relmajor = `pdv_relmajor_expected_mean_avg'/((1-$r)^`t')
						  	local pih_quant_major = `pdv_major_expected_quant_avg'/((1-$r)^`t')
						  	local pih_quant_relmajor = `pdv_relmajor_expected_quant_avg'/((1-$r)^`t')
							
							
							matrix  pih = pih\(`t' + 20, `pih_mean_major',`pih_quant_major',`pih_mean_relmajor',`pih_quant_relmajor')
					   }

						matrix pih = pih[2...,1...]
				   
				   *CALCULATE UTILITY - MEAN
					 local utility_major = 0
					 local utility_relmajor = 0

					 forvalues t = 1/$numquarters {
						if pih[`t',2] > 0  local utility_major = `utility_major'+ `beta'^(`t')*(pih[`t',2]^(1-`gamma') )/(1-`gamma')
						if pih[`t',4] > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(pih[`t',4]^(1-`gamma') )/(1-`gamma')
					 }
					 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
					 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (0) ("Permanent Income Hypothesis") (`utility_major') (`utility_relmajor') (`rel_utility')
	
				   *CALCULATE UTILITY - QUANTILES
					 local utility_major = 0
					 local utility_relmajor = 0

					 forvalues t = 1/$numquarters {
						if pih[`t',3] > 0  local utility_major = `utility_major'+ `beta'^(`t')*(pih[`t',3]^(1-`gamma') )/(1-`gamma')
						if pih[`t',5] > 0  local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(pih[`t',5]^(1-`gamma') )/(1-`gamma')
					 }
					 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
					 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (1) ("Permanent Income Hypothesis") (`utility_major') (`utility_relmajor') (`rel_utility')
				   

		
				   
		*SAVINGS RULE - SAVE ANYTHING ABOVE KNOWN INCOME STREAM, DRAW DOWN SAVINGS UNTIL HIT PREDICTED INCOME OR SAVINGS = 0
		
			 *MEAN
			 local utility_major = 0
			 local utility_relmajor = 0
			 local savings_major = 0
			 local savings_relmajor = 0
			 forvalues t = 1/$numquarters {
				if `savings_major' >= 0 { 
					 local consumption_major = min(predicted_earnings[`t',2], predicted_earnings[`t',3] + `savings_major')
				}
				else {
					 local consumption_major = predicted_earnings[`t',3]
				}
			    local savings_major = (predicted_earnings[`t',3] + `savings_major' - `consumption_major')*(1+$r) 
				
				if `savings_relmajor' >= 0 { 
					 local consumption_relmajor = min(predicted_earnings[`t',7], predicted_earnings[`t',8] + `savings_major')
				}
				else {
					 local consumption_relmajor = predicted_earnings[`t',8]
				}
			    local savings_relmajor = (predicted_earnings[`t',8] + `savings_relmajor' - `consumption_relmajor')*(1+$r) 

				
				if `consumption_major' > 0 local utility_major = `utility_major'+ `beta'^(`t')*(`consumption_major'^(1-`gamma') )/(1-`gamma')
				if `consumption_relmajor' > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(`consumption_relmajor'^(1-`gamma') )/(1-`gamma')
			}
				
			 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
			 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (0) ("Savings Rule") (`utility_major') (`utility_relmajor') (`rel_utility')

			 *QUANTILE
			 local utility_major = 0
			 local utility_relmajor = 0
			 local savings_major = 0
			 local savings_relmajor = 0
			 
			  forvalues t = 1/$numquarters {
				if `savings_major' >= 0 { 
					 local consumption_major = min(predicted_earnings[`t',4], predicted_earnings[`t',5] + `savings_major')
				}
				else {
					 local consumption_major = predicted_earnings[`t',5]
				}
			    local savings_major = (predicted_earnings[`t',5] + `savings_major' - `consumption_major')*(1+$r) 
				
				if `savings_relmajor' >= 0 { 
					 local consumption_relmajor = min(predicted_earnings[`t',9], predicted_earnings[`t',10] + `savings_major')
				}
				else {
					 local consumption_relmajor = predicted_earnings[`t',10]
				}
			    local savings_relmajor = (predicted_earnings[`t',10] + `savings_relmajor' - `consumption_relmajor')*(1+$r) 
	
				if `consumption_major' > 0 local utility_major = `utility_major'+ `beta'^(`t')*(`consumption_major'^(1-`gamma') )/(1-`gamma')
				if `consumption_relmajor' > 0 local utility_relmajor = `utility_relmajor'+ `beta'^(`t')*(`consumption_relmajor'^(1-`gamma') )/(1-`gamma')
			}
				
			 local rel_utility = (`utility_major'- `utility_relmajor')/abs(`utility_relmajor')			
			 post simulation_2yr ("`major'") ("2yr") (`beta') (`run') (`gamma') (1) ("Savings Rule") (`utility_major') (`utility_relmajor') (`rel_utility')

		
	*CLOSE LOOP FOR RUN
	}
*CLOSE LOOP FOR MAJOR
}

*CLOSE LOOP FOR CRRA (GAMMA)
} 	

*CLOSE LOOP FOR DISCOUNT RATE (BETA)
}	
		postclose simulation_2yr
		use "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_2yr.dta" , clear
		save  "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_2yr_runs.dta", replace


		
		*COLLAPSE TO SUMMARY DATA BY MAJOR AND TYPE
		collapse (median) utility_major utility_libarts relutility, by(major scenario gamma quantiles discountrate)
		sort major scenario quantiles  gamma 
		matrix drop _all
		save  "C:\Users\imberman\Dropbox\Texas Majors\Major Earnings Trajectories and Variability\Results\Utility model\simulation_results_quantiles_2yr.dta", replace





